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Abstract 

We have made a study of several update algorithms using the XY 
model. We find that sequential local overrelaxation is not ergodic 
at the scale of typical Monte Carlo simulation time. We have intro- 
duced a new multisite microcanonical update method, which yields 
results compatible with those of random overrelaxation and the mi- 
crocanonical demon algorithm, which are very much slower, all being 
incompatible with the sequential overrelaxation results. 
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Microcanonical local algorithms are used often in numerical simulations of 
spin systems or lattice gauge theories. Sometimes they are used in combina- 
tion with canonical algorithms in order to accelerate decorrelation in Monte 
Carlo time. Some other times one is interested in a pure microcanonical 
simulation for physical reasons |jl|, |^. We do not consider nonlocal mi- 
crocanonical methods which are based on molecular dynamics and are used 
specially for studying dynamical fermionslQ. 

Possibly, the most popular microcanonical algorithm is Local Microcanon- 
ical Overrelaxation (LMO) The algorithm, simple and fast, has a dynam- 
ical exponent z = 1 when the update is done sequentially, and z = 2 when it 
is done randomly Q], the difference being due to a wave effect occurring in 
the former case, which causes the change made in the update of a variable to 
propagate to the next one, the updated configuration being farther from the 
original one than in the case consecutive updates had been made indepen- 
dently. Yet, no proof exists of the algorithm being ergodic. Its determinism, 
particularly in the case of sequential update, might induce one to think that 
it is not ergodic. Indeed, we shall present evidence of the lack of ergodicity 
of the sequential LMO method. 

We shall introduce here the Multisite Microcanonical method (MM), a 
new microcanonical updating algorithm which exchanges energy among dif- 
ferent regions of the lattice at each update. It has two features which are 
absent from the LMO: it is non-local and non-deterministic. They contribute 
to diminishing Monte Carlo correlation time, and avoid the biases that plague 
LMO. 

Apart from the mentioned LMO and MM, we have also done simulations 
using the Microcanonical Demon Algorithm, (MDA), introduced by Creutz 
0, in order to have one further point of reference. 

We have chosen a simple model for our simulations, the two dimensional 
XY model. We remark that the method could be easily generalized to 0{N) 
models in arbitrary dimensions, or even to abelian and SU(2) gauge theories. 

In what follows, we shall recall the two dimensional XY model, introduc- 
ing the observables we have measured, and study the mentioned algorithms, 
LMO, MDA and MM. 

In a square lattice of side L, with periodic boundary conditions, we define 
on each site a variable G U{1), the action 
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and the density of states with energy E 

N{E) = [ lld(l>{n)6{S{<l>) - E). (2) 

■' n 

Z being the partition function. 

The presented results have been obtained in a 64^ lattice, with a nor- 
mahzed energy value E — 0.4453, which corresponds to /3 ~ 1.0. The data 
plotted in figures have been obtained with 45, 000 configurations (after ther- 
malization). Where the results were not conclusive we have added 500,000 
MC iterations. Most runs have been started from two fixed (thermalized) 
configurations. Errors have been calculated using the jackknife method. 

Among the measured observables, the energy, the process being micro- 
canonical, should remain constant. We have measured it as a check of the 
different methods, observing small changes due to rounding errors in the mi- 
crocanonical methods, and larger oscillations in the MDA case, because in 
the latter only the total energy (system plus demon) remains constant. 

We have measured the total magnetization, defined as M = -p- Y,n 
M being a complex number, we have measured its modulus and phase, which, 
although it can be chosen arbitrarily by a global gauge change, 0(n) U(j){n), 
Vn, U e U(l), can give us a hint on how the system is evolving. 

The most useful observables for our study have turned out to be the 
spatial correlations. In particular, we have measured the correlation between 
parallel lines at distances from to L/2, i.e., defining 

^(^) ^ T 5Z ^(^' y) ^(y) = 7 '^(^' y) (3) 

^ y ^ X 

(where the notation (f){x,y) has been used, instead of 4>{n)), the correlation 
at a distance r is written 

rj{r) = ^{uj{x)uj*{x + r)+ u{y)u*{y + r)). (4) 

We have studied separately the real and imaginary parts of r]{r). From the 
real part, correlation lengths can be extracted. 

The imaginary part should be zero at any r, due to the symmetry of 
the action under the transformation 0(n) — > 0*(n),Vn. We shall see that 
the imaginary part is not zero in the LMO method. This implies that the 
system does not evolve uniformly in the whole phase-space, so breaking the 
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fundamental assumption of the microcanonical ensemble: at any time, all 
the states of the system with the given energy have equal probability. 

Let us consider first the LMO algorithm. Writing = e*^'-"^ and defining 

pK)e*"^"«) = 5]e-*'("°+'^), (5) 

the action associated to the site no reads 

S{no) = Re (0(no) J2 <^*('^o + At)) = p(^o) cos(6'(no) + a('^o))- (6) 

One LMO step consists on picking a site hq and updating 6{nQ) to a new 
value 

9'ino) = -9ino) - 2a(no) (7) 

which, due to the parity of the cos function, does not alter the energy. 

The step satisfies detailed balance in a deterministic way: the probability 
of going from one state to the other is 1 in both senses. But a simulation 
takes many of these steps, and the order in which they are taken is also 
important. If we question the relationship between two states separated by N 
such steps, taking the intermediate sites randomly satisfies detailed balance, 
but taking them sequentially does not. The looser condition of balance |^ 
is satisfied, though. The latter, together with ergodicity, is sufficient for the 
algorithm to reproduce the desired distribution after a reasonable number of 
sweeps. However, the determinacy of each single step makes one question 
the ability of the method to satisfy ergodicity. Running a sweep forward 
and backwards would leave the configuration unchanged. Also, a sequential 
overrelaxation sweep on a one- dimensional XY model would have as only 
effect to transport a little energy packet along the chain, leaving the rest of 
the chain undisturbed. 

Proving ergodicity is, in general, difficult, but disproving it may be eas- 
ier. A proof of lack of ergodicity is the occurrence of results which, in the 
limit of infinite statistics, depend on the initial state. Also, finding a con- 
served quantity which should not be conserved means that the algorithm 
runs in a certain region of phase space, characterized by a certain value of 
the conserved quantity, without being able to leave it. 

Let us consider sequential LMO, SLMO. Fig. la shows that the phase 
runs over the interval (— vr, tt] at a practically constant rate (we shall call 
that property helicity). The sign and modulus of the helicity depend on the 
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Figure 1: Evolution in Monte Carlo time of the phase of the magnetization 
for several algorithms. Only a sample of each run is shown. 
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initial configuration, and they are kept through the simulation with small 
fluctuations. Wc have consistently used two fixed configurations, character- 
ized by their respective large and small value of the helicity, as initial con- 
figurations in all our simulations, in order to check the possible dependence 
of the results on the initial state. Making the necessary allowances for the 
proviso made above on the meaning of the phase, the origin of the rotation 
described being a global rotation of the system would be of no importance, 
but problems might arise if the rotations were local. 

The imaginary part of the spatial correlation defined above, Im(r;(r)), 
should average statistically to zero, the quantity being invariant under a 
global rotation, but sensitive to local ones. Fig. 2a shows that SLMO yields 
values for Im(?7(r)) clearly incompatible with zero within the statistical er- 
rors, statistics being sufliciently high. That, together with the fact that other 
algorithms do yield results fully compatible with zero, confirms our suspicion 
that SLMO does not run over the phase space properly. 

As Im(77(r)) is not either a very usual observable, we have considered too 
Re(?7(r)). Fig. 3a shows that SLMO yields results incompatible with each 
other and dependent on the helicity of the initial configuration. In order to 
reinforce this conclusion we have increased the statistics with 500, 000 more 
iterations, with results confirming the incompatibility. 

Randomly ordered LMO (RLMO) does not conserve helicity. In that case, 
the phase of M, shown in Fig. lb, evolves very slowly in a nonsystematic way. 
Fig. 2b shows that lm{r]{r)) is compatible with zero. However, Fig. 3b shows 
that the real part is also dependent on the initial configuration, although 
the values are less different from each other than in the SLMO case, and 
the errors, also for |M|, are much bigger, and so are the correlation times. 
Clearly, it would be wrong to infer from the error size that SLMO is better 
than RLMO, as we know by now that the former runs only over a small 
region of phase space. 

Next, we have run also up to 500, 000 iterations for each initial configu- 
ration (high and low helicity) and in this case the results for Re{r}{r)) have 
changed dramatically with respect to those obtained with 50, 000 iterations. 
Now the results for the two different starts become fully compatible (and 
compatible also with those obtained by other methods shown later). There- 
fore RLMO behaves ergodically, but with a Monte Carlo time scale very large 
compared with the scale of the other algorithms shown later. 

It must be said, in discharge of SLMO, that despite its lack of ergodicity, 
the fact that it decorrelates efficiently, in the sense that it produces configu- 
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Figure 2: Imaginary part of the spatial correlation evaluated by several al- 
gorithms (with 50, 000 iterations in each case). 



rations which, although all belong to a restricted region of phase space, are 
fairly far from each other, has been widely shown, make it very use- 

ful when combined with other algorithms, like Metropolis, which reinstate 
ergodicity. 

Let us consider next the Microcanonical Demon Algorithm, in which an 
extra degree of freedom, the demon, travels through the system, transferring 
energy from one point to another. The total energy of the system plus 
the demon is kept constant, which makes the algorithm to be not strictly 
microcanonical from the point of view of the system alone. The algorithm 
works as follows: a small amount of absolute energy is assigned initially to 
the demon, say, 1. When updating a spin, a new random tentative value 
is chosen, and the new system energy is computed. In case it increases, 
the change is accepted, and the excess energy is given to the demon. If the 
energy decreases, the change is accepted only if the demon energy is sufficient 
to provide the difference, the demon energy being always kept non-negative. 

As only one demon is being used (more demons can be included in the 
formalism) and its initial energy is small, the system evolves in a practically 
microcanonical way. Fig. Ic shows the evolution of the phase of M, which 
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Figure 3: Real part of the spatial correlation evaluated by several algorithms 
(with 50, 000 iterations in each case). 



shows no helicity. Fig. 2a shows lm{ri{r)) calculated by the MDA. As can be 
seen, it is compatible with 0. In Fig. 3c we plot the results of the algorithm for 
Re(77(r)) for the two fixed initial configurations, showing the independence 
of the evolution from the initial state. 

Next, we introduce our Multisite Microcanonical method, which is a non- 
local, non-deterministic algorithm allowing the exchange of energy among 
different sites, keeping the system energy constant. 

Using the notation introduced in we rewrite the system energy as 

E = lReY, 0(n)0*(n + /x) = Wpi^) cos(^(n) + a(n)) = ^(p,f) . (8) 

Keeping a constant energy means keeping constant the scalar product of 
the two vectors p and x just defined: 

(p, i) = K (9) 

The expression is valid for the energy of a subset of any iV sites, the 
vectors p and x being now iV-dimensional. If, moreover, the set does not 
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contain interacting sites, updating the set means changing the vector x, as p 
and the a's remain constant in the process. 

The algorithm proceeds as foUows: N non-interacting sites are chosen, p 
and X are computed, a new x' is generated in such a way that (p, x') = (p, x) 
and the new variables 9' are obtained from x'. 

Yet, many technical details must be kept in mind. The components of the 
vector x' being cosines, must lie in the interval [—1,1], and, when computing 
the new values 6', one has the freedom to choose the sign of sin(^'-|-Q;), main- 
taining detailed balance. One possibility is choosing it with equal probability 
for both signs, another is to include a measure of overrelaxation, choosing 
the sign which yields the 9' which lies farther away from the original one. 

But, in order to ensure ergodicity and balance, it is the 0' space (or 
equivalently, the 6' space) which must be filled uniformly. Once we have 
generated the x's uniformly, we can get a uniform distribution in the 9' space 
weighting the configurations with the Jacobian of the transformation from 
the ^'s to the x's, which is 

J{9 -^x) = r\ . ^ (10) 

^ ^ V sm(^(n) + Q;(n)) ^ ^ 

so that, once we have generated uniformly the new x' we keep it or not with 
a probability proportional to the quotient of the new and old jacobians. 

Our problem reduces, then, to uniformly generating points in a (N — 1)- 
dimensional polyhedron, which is the intersection of the A^- dimensional solid 
hypercube of side two centred in the origin and the (A^ — l)-dimensional 
hyperplane perpendicular to p containing x. 

An obvious way to generate x uniformly is generating uniformly {N — 1) 
components in the interval [1,-1], and finding the A^^th component which 
ensures that x lies in the hyperplane. In case the last component lies in the 
interval [1,-1], one accepts the new point, otherwise one must start and try 
again. Clearly, the inefficiency of this method grows exponentially with N, 
making the method useless for other than a few sites. 

The largest set of non-interacting sites is the set of all equally coloured 
sites in the checkerboard lattice, which contains N = V/2 points, which for 
an interesting lattice is of the order of thousands, which makes the use of 
such inefficient methods hopeless. 

We have explored other more sofisticated methods, on which we shall 
report elsewhere, but in general one must trade acceptance for computational 
load, the net result being the unimplementability of the method for large N. 
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Figure 4: Intersection of the solid hypercube of side 2 centered in the origin 
with the hyperplane perpendicular to p containing x in two dimensions. 



We have relaxed the demand to generate uniformly at each step, still insisting 
on keeping detailed balance, but then the explored region becomes small, and 
the decorrelation is poor. 

Turning to MM with small N, N — 1 \s local and corresponds to LMO. 
N — 2\s the smallest value which allows a non-local update, and consequently 
an energy exchange between arbitrarily different sites (LMO exchanges en- 
ergy only between the ends of a link). That, together with the fact that at 
N = 2 the efficiency of the method is highest, the geometrical acceptance 
rate being one, as will be shown later, has moved us to choose N — 2. 

Having chosen the value of N, the problem remains of how to pick the two 
non-interacting sites. The most unbiased choice is picking them at random 
(RMM), and that has been done. Yet, we have mentioned that sequential 
updating has the advantage of a wave effect, which transmits the decorrela- 
tion effect along the sequence, so we have basically done simulations moving 
one site sequentially and picking the other one at random, which seems to 
be efficient. 

With N — 2, the equation 

pixi + P2X2 = K (11) 
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must be satisfied, Xi and X2 being uniformly distributed in tlie allowed region. 
That is easily done; one looks for the maximum and minimum values of Xi 
for which a solution for X2 to (|Tl]) exists between —1 and 1, and then a value 
for xi is generated randomly and uniformly in that interval. Fig. 4 illustrates 
the process. 

Let us call {xi,X2) the old values of x and {x'i,x'2) the new ones, gen- 
erated as described above. From them the values of the new {d'1,9'2) are 
obtained, solving x' = cos{6' + a). Next, in order to obtain uniformly dis- 
tributed ^'s, the quotient C of jacobians is computed 

^ ^ {sm{e\ + a\)sm{e'2 + a'2)r' ^^3) 
(sin(^i + ai) sin(^2 + "2))" 

and the new values are accepted if C > 1, or otherwise with probability 
proportional to C. The acceptance from this factor becomes around 70%. 
Finally we have chosen to take the sign of sin(^' + a) which yields the farthest 
6' from the original one. 

The results of the simulation with the implementation of the MM algo- 
rithm just described, are as follows: Fig. Id shows the evolution of the magne- 
tization phase, which rambles unsystematically. Fig. 2b shows that lm{ri{r)) 
is compatible with 0, with much smaller errors than for other algorithms. 
Fig. 3d , in which results of RMM runs are included, shows that Re(r7(r)) 
behaves identically for both initial configurations and both implementations 
of the algorithm, with smaller errors for the choice sequential-random. 

MM results for Iie{7]{r)) are just outside error bars with respect to MDA 
results. This is due to the fact, already mentioned, that MDA simulations 
are not strictly microcanonical, and in fact the energy values of the MDA 
simulation, which run in the range (0.4550,0.4562), yield a slightly higher 
mean value than those of the MM simulations, which range in (0.4552, 0.4553) 
(due to rounding). This explains qualitatively why the MDA results are 
slightly bigger than the MM ones. This situation is unavoidable, if we insist 
on studying the dependence of the results on the initial state, which forces 
us to use fixed initial states for all our simulations, so losing in a certain 
measure our control on the average energy in the case of MDA. 

Forgetting about helicity, we have run with MM at energy 0.4556 (the 
mean energy of MDA) and in this case the results are compatible. 

The comparison with the LMO results in Fig. 3a,b is pointless, since those 
results are selfincompatible. 
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Our conclusions are, then: 

SLMO is not ergodic (for the two-dimensional XY model). 

We have introduced a new multisite, microcanonical update method, MM. 

Among the microcanonical algorithms studied here, SLMO and MM have 
proved to be the most efficient, as far as Monte Carlo time decorrelation is 
concerned, although LMO runs only over a restricted region of phase space. 

The RLMO, MDA and MM algorithms yield compatible results, while 
LMO results are start-dependent and incompatible. 

We wish to thank Juan J. Ruiz-Lorenzo for useful discussions. Partially 
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